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Abstract. We compute polarization maps for molecular cloud cores modeled as 




magnetized singular isothermal toroids, under the assumption that the emitting 




dust grains are aspherical and aligned with the large-scale magnetic field. We show 




that, depending on the inclination of the toroid with the line-of-sight, the bend- 




ing of the magnetic field lines resulting from the need to counteract the inward 




pull of gravity naturally produces a depolarization effect toward the centre of the 




map. We compute the decrease of polarization degree with increasing intensity for 




different viewing angles and frequencies, and we show that an outward increasing 




temperature gradient, as expected in starless cores heated by the external radi- 




ation field, enhances the decrease of polarization. We compare our results with 




recent observations, and we conclude that this geometrical effect, together with 




other mechanisms of depolarization, may significantly contribute to the decrease 




of polarization degrees with intensity observed in the majority of molecular cloud 




cores. Finally, we consider the dependence of the polarization degree on the dust 




temperature gradient predicted for externally heated clouds, and we briefly com- 




ment on the limits of the Chandrasekhar-Fermi formula to estimate the magnetic 




field strength in molecular cloud cores. 
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. Introduction 





Mapping the polarization of the thermal emission of dust at millimitre or submillime- 
tre wavelengths (usually A = 850 /im or 1.3 mm) is the principal means of probing 
Send offprint requests to: J. Gongalves, 
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the magnetic field geometry in molecular cloud cores. A frequent characteristic of these 
observations is the decrease of polarization degree p as a function of the total observed in- 
tensity /. Usually p decreases with increasing intensity / with a power-law behavior, from 
a maximum value of ~ 15% to about the observable limit of ~ 1%. This depolarization 
effect (sometimes referred to as "polarization hole" , or "polarization limb brightening" ) 
has been observed in many dense cores and filamentary clouds. Recent examples in- 
clude the OMC-3 region of the Orion A filamentary molecular cloud where p cx i~ - 7 
(Matthews & Wilson 2000), several dense cores in the dark cloud Barnard 1, where 
p oc J -0 8 (Matthews & Wilson 2002), Bok globules mapped by Henning et al. (2001) 
where p oc I~~ o e , and dense cores mapped interferometrically by Lai et al. (2002), where 
p oc I~ o s . In some cases, the decrease of polarization with intensity is quite steep: for 
example, Crutcher et al. (2004) find p oc in the dark cloud L183, implying that not 

only the polarization degree but also the polarized intensity I p — pi decreases toward 
the centre of this cloud. A similar situation is apparently found also in the starless core 
L1544 (Ward-Thompson et al. 2000). 

The dark cloud Barnard 1 is an excellent example of the type of object of interest 
to us. This cloud, mapped with SCUBA by Matthews & Wilson (2002), has dimensions 
of roughly 0.2 x 0.4 pc, densities of a few times 10 4 cm -3 , and temperatures of order 
10 K. These appear unfavorable characteristics for aligning grains (see e.g. Lazarian et 
al. 1997). Nevertheless, Bl exhibits extremely uniform polarization at a level of around 
5% and with an extremely uniform position angle. This last point in our opinion is a 
strong argument in favor of an essentially uniform and likely dominant magnetic field 
over most of the core. However, there are exceptions to the rule in the shape of four 
high density (of order 10 6 cm~ 3 ) prestellar or protostellar "inclusions" of dimensions 
~ 0.03 pc (Bl-a,b,c,d) apparently embedded within the region of the "general core" 
referred to above. These extremely high density cores also surprisingly show evidence for 
polarization, albeit at a lower level and with different position angles than the general 
core. In particular, in the core Bl-c the polarization degree is very low and almost uniform 
over the whole observed intensity range (Matthews & Wilson 2002) 1 . For illustration 
purposes, the p-I relations measured in three well studied cores (L183, L1544, and Bl-c) 
are collected in Fig. ^ together with power-law fits of the original data. 

The observed decrease of polarization with increasing intensity has been attributed 
to changes in the dust grain optical properties or shapes in the cold cores interiors due 
to grain-growth (Vrba et al. 1993, or Hildebrand et al. 1999): not only are bigger grains 

1 We note however parenthetically that the observations of polarization at high densities (n ~ 
10 6 cm -3 ) where radiation fields can hardly be of importance represent a serious challenge to 
grain alignment theorists. In fact, recent work shows that radiative torques due to starlight are 
required to drive grains to the suprathermal rotation rates necessary to minimize the disaligning 
effect of random collisions (Draine & Weingartner 1996, 1997). 
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Fig. 1. Distribution of percentage of polarization p vs. intensity I/I max at A = 850 /xm, 
normalized to its maximum value, for the starless cores L183 (filled triangles, from 
Crutcher et al. 2004), L1544 (empty triangles, from Ward-Thompson et al. 2000), and 
the Bl-c core (circles, from Matthews & Wilson 2002). The dashed lines are power-law 
fits, with slopes -1.6 for L1544, -1.2 for L183, and -0.4 for Bl-c. 



more difficult to align than smaller grains, but also the agglomeration process may make 
grains more spherical, thus further reducing the alignment efficiency. 

All these effects are difficult to quantify. A larger density tends in general to dealign 
grains (Lazarian et al. 1997), but the enhanced magnetic field strength should have 
the opposite effect. On the basis of the Davis & Greenstein (1951) alignment model by 
paramagnetic relaxation, Vrba et al. (1981) showed that the ratio pi Ay (equivalent to 
p/I at millimetre wavelengths) should scale as B 2 a~ 1 n~ 1 , where B is the intensity of the 
magnetic field, a the grain size, and n the ambient density. If the magnetic field strength 
scales as n 1 ! 2 in molecular clouds, as empirically observed (Crutcher 1999), the effect of 
the increased density and magnetic field cancel, and one is left with the dependence of 
pjAy on the grain's size (and shape). 

Fiege & Pudritz (2000) showed that a helical magnetic field geometry in a cylindri- 
cal cloud with uniform grain properties could naturally produce a depolarization effect, 
in agreement with the submm polarization observations of the OMC-3 region of the 
Orion A filamentary molecular cloud (Matthews & Wilson 2000; Matthews, Wilson, & 
Fiege 2001). Padoan et al. (2001) were able to produce a decreasing p vs. I relation 
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in three-dimensional MHD simulations of supersonic and super- Alfvcnic turbulence as- 
suming that grains are not aligned above Ay ~ 3 mag, but the dynamic range of the 
simulation did not extend much beyond this value of extinction. They noticed, however, 
that even assuming a uniform grain alignment efficiency, a decrease of p with / could be 
reproduced for particular orientations of the core magnetic field relative to the line of 
sight, an effect that will be further explored in the present paper. 

It is possible in principle that the observed depolarization may be accounted for 
(at least in part) by beam smearing over tangled, small-scale field structures: Rao et 
al. (1998), for example, found that the decrease in polarization toward the Klcinman-Low 
nebula previously seen with single-dish observations was a result of subrcsolution-scale 
variations in the magnetic field that are averaged out by larger beams. However, given 
the relatively high levels of polarization detected at the core peaks, this effect should not 
be dominant. 

We will try in this paper to further analyze the influence of the magnetic field geometry 
on the observed depolarization effect, studying cloud models dominated by a large scale 
magnetic field. In particular, we show that the pinching of the magnetic field expected 
in dense, self-gravitating molecular cloud cores naturally produces a decrease of the 
polarization degree toward the centre of the core, for a large range of viewing angles. 
In addition, we find that this geometrical depolarization effect is further enhanced by a 
dust temperature gradient increasing outward, as expected in externally heated starless 
cores (Evans et al. 2001, Zucconi et al. 2001; Stamatellos & Whitworth 2003; Gongalves 
et al. 2004, hereafter GGW). However, while these effects can contribute substantially to 
the observed depolarization, the measured p—I relations are sometimes steeper than our 
predictions. Thus, we do not suggest that the field geometry and the dust temperature 
distribution are the sole means by which these observations could be explained: variations 
in the grain optical and/or geometrical polarization properties in the densest parts of a 
cloud are still required to account for the full range of observed polarization values. 

The structure of the paper is the following: in Sect. 2, we describe the model adopted 
for magnetized molecular cloud cores and define the relevant physical and geometrical pa- 
rameters of the problem; in Sect. 3, we show synthetic polarization maps computed with 
these models, and illustrate the effect of geometrical depolarization as well as its depen- 
dence on inclination, wavelength, and temperature distribution. In Sect. 4 we summarize 
our conclusions. 

2. The model 

We model molecular cloud cores as singular isothermal toroids, i.e. scale-free, axisym- 
metric equilibrium configurations of an isothermal gas cloud under the influence of self- 
gravity, gas pressure and magnetic forces (Li & Shu 1996). These toroids are characterized 
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by one non-dimensional parameter, Ho, representing the fractional amount of support 
provided by magnetic forces. We show our polarization results for three models, with 
H = 0.2, Hq = 0.5 and H = 1.25, corresponding to mass-to-flux ratios (in units of the 
critical value) A ks 6, 3, and 1.6. With increasing Ho, the density distribution becomes 
flatter, and the configuration becomes a thin disk for Hq — > oo. Figure [2] shows isodensity 
contours and magnetic field lines for a singular isothermal toroid with Hq = 1.25, and 
illustrates some geometrical quantities adopted in our analysis. 




Fig. 2. The geometry of the model. The Hq = 1.25 singular isothermal toroid (thick 
curves: isodensity contours; thin curves: magnetic field lines) is observed from a line of 
sight (l.o.s.) inclined by an angle 9 with respect to the equatorial plane of the toroid 
(9 = 0: edge-on; 6 = 90°: pole-on). Dust grains are assumed to be aligned with the local 
direction of the magnetic field, making an angle 7 with the plane of the sky (p.o.s.). 
Notice that for the inclination shown in this figure, magnetic field lines are almost in the 
plane of the sky for lines of sight intercepting the toroid above and below the centre, 
resulting in a relatively higher polarization degree of the outer parts of the core with 
respect to the central region. 

We further assume that the core is bathed in the average interstellar radiation field 
(ISRF) of the solar neighborhood (Black 1994), and use the results of GGW for the 
temperature distribution. The dust temperature profile and the intensity of the emitted 
radiation at submillimetre wavelengths for the Hq — 0.5 toroid were presented in Sect. 5 
and Fig. 3 of GGW. We note some inconsistencies in our approach: first, the core's 
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density distribution was obtained by Li & Shu (1996) under the hypothesis that the gas 
is isothermal, whereas we explicitly consider deviations from isothermality at least in the 
dust component; second, GGW computed the dust temperature distribution assuming 
spherical grains, while in the present work, for the purpose of computing the polarization 
of the emitted radiation, we assume that the grains are aspherical and aligned with 
the magnetic field. However, the results of GGW arc not substantially affected by grain 
shape. 

Since we are interested in computing the sub-millimetre thermal emission, we neglect 
self-absorption of radiation. We also assume that the gas-to-dust ratio and the proper- 
ties of the dust grains with respect to absorption and polarizing efficiency are uniform 
throughout the cloud. The polarization degree p and polarization angle \ (the direc- 
tion of polarization in the plane of the sky) are defined in terms of the standard Stokes 
parameters Q, U and /, 



tan2 X =^. (2) 

Here, we compute p and \ following a method developed by Lee & Draine (1985), and 
elaborated by Wardle & Konigl (1990), Fiege & Pudritz (2000), and Padoan et al. (2001), 
but we generalize their method allowing a dependence of the dust temperature upon 
position inside the core. In this formulation, p and x are given by the simple expressions 

tan2 X =-, (4) 
<7 

where a is a nondimensional parameter representing the "polarizing efficiency" of dust 
grains, q and u are the "reduced" Stokes parameters 

q = j pB x (T d ) cos 2tp cos 2 7 dl , (5) 

u = J p J B A (T d )sin2V>cos 2 7 dl, (6) 
and 



S = J P B x {T d ) dl, (7) 

S 2 = j pB x (T d ) - I) dl. (8) 

These four quantities are integrals along the line of sight I of the product of the total 
density p (proportional, under our hypothesis, to the number density of dust grains), the 
Planck function B\(T d ) (representing the dust emissivity at the dust temperature T d ), 
and a geometric factor accounting for the orientation of the magnetic field at each point 
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inside the cloud, characterized by the two angles ip and 7. Specifically, the former is the 
angle between a direction in the plane of the sky (e.g. north) and the component of B in 
that plane; the latter is the angle between the local direction of B and the plane of the 
sky (see Fig. 0. 

Thus, knowledge of the density and temperature distributions p(t) and Td(£), together 
with the geometry of the magnetic field in a cloud core expressed by i/)(£) and 
are sufficient to completely determine the polarization characteristics of the radiation 
emitted by dust grains, if one specifies the nondimensional parameter a (constant under 
our assumptions), containing all information about the absorptions cross sections and 
the alignment efficiency. The numerical value of a can be easily fixed observing that the 
maximum polarization degree is achieved when B is in the plane of sky. In this case, with 
V> constant, eq. 10 gives (Fiege & Pudritz 2000) 

6a . . 

Pmax = 7 , (") 

— a. 

or 

6p m ax . c 1 f-\t\\ 

a = — m p m ax II Pmax < 1- (10) 

u t Pmax 

The choice of a is therefore equivalent to a normalization of the expected degree of po- 
larization. The polarization degree and angle can then be computed fixing the maximum 
polarization degree, instead of assuming specific grain properties. In this work we choose 
a = 0.15, which corresponds, from eq. 10, to p max ~ 15%. 

3. Results 

For a molecular cloud core modeled as described in Sect. 2, with specified density and 
temperature distributions, the only free parameters are a and the inclination angle 8 
of the line of sight with respect to the equatorial plane of the toroid (8 — 0: edge-on; 
8 = 90°: pole-on). In this section we show the polarization of the emitted dust radiation 
for a = 0.15, varying the inclination of the toroid (8 — 20°, 40°, and 60°), and for toroids 
with H = 0.2, 0.5 and 1.25. 



3.1. Geometric depolarization at intermediate inclination 

In Fig. we show the isophotes and polarization vectors (rotated by 90°) of the dust 
emission computed for the H = 0.2 toroid at A = 850 fxm for = 20°, 40° and 60°. The 
model results have been convolved with a telescope beam of FWHS of 12", assuming a 
distance of 150 pc in all cases. Figures 0] and [S] show the same maps for the H = 0.5 and 
1.25 models, respectively. The polarization vectors have been rotated by 90° to show the 
approximate average orientation of the magnetic field in the plane of the sky. 

The polarization pattern is clearly symmetric, because of the assumed axial symme- 
try of the model. The equatorial pinching of the magnetic field lines (see Fig. |2J has 
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Fig. 3. Maps of the dust emission and polarization at A = 850 ^tm for the Hq = 0.2 
singular isothermal toroid, shown in the bottom right panel (solid curves, isodensity 
contours, dashed curves, magnetic field lines). The first three panels are for inclination 
with respect to the plane of the sky 6 — 20° (top left panel), 9 = 40° (top right panel), 
and 9 — 60° (bottom left panel). Each vector is proportional to the polarization degree 
(see scale on the lower right corner of each panel), and has been rotated by 90° to show 
the average orientation of the magnetic field in the plane of the sky. The intensity is shown 
by contours logarithmically spaced by 0.2 dex starting from 10% of the peak value. 

an important effect on the non-uniformity of the polarization degree across the cloud, 
resulting in a significant decrease of p toward the central regions. In fact, for a hourglass 
magnetic field configuration observed at intermediate inclinations (9 « 30°- 40°), the 
largest component of the field in the plane of the sky is found in the outer parts of the 
cloud, whereas lines of sight close to the cloud's centre intercept regions where the bend- 
ing of field lines is stronger and therefore the component of the field in the plane of the 
sky is relatively weaker (cf. Fig. EI> ■ 

The depolarization effect is better illustrated in Fig. H3 showing the degree of polar- 
ization p as function of the intensity at 850 /im (normalized to the peak intensity) for the 
H = 0.2, 0.5 and 1.25 toroids at three inclinations, 9 = 20°, 40° and 60°. The minimum 
level of polarization is usually attained toward the centre of the core, where the non 
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Fig. 4. Same as Figure|2|for the Ho = 0.5 singular isothermal toroid. 

uniformity of the magnetic field along the line of sight is larger, and therefore cancel- 
lation effects more important. Over about one order of magnitude increase in intensity, 
the polarization degree decreases with roughly a power-law behavior. It is interesting to 
notice that even for a moderate pinching of the field (e.g. for the Ho = 0.2 toroid), the 
decrease of polarization towards the centre is already significant. 

We do not attempt in this paper to model specific objects, but we stress the qualitative 
similarity between our theoretical p-I diagrams shown in Fig. and the p-I relations 
observed in dense cores (Fig. In some cases, however, like in the dark cloud L183 
(Crutcher et al. 2004) or L1544 (Ward-Thomson et al. 2000), the observed p-I relation 
is too steep to be explained only on the basis of the field morphology and inclination. 
A full model for the polarization of the emitted radiation probably requires additional 
mechanisms (such as an increase in size and sphericity of dust grains near the core centre) 
to reduce the polarizing efficiency a at high values of density or extinction. 

3.2. Dependence of depolarization on aspect ratio 

The maximum and minimum value of the polarization degree, p max and p m i n , respectively, 
depend on the inclination angle 9 of the core with respect to the line of sight. Both p max 
and Pmin decrease when the inclination of the toroid changes from edge-on (6 = 0) to 
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Fig. 5. Same as Figure and 0] for the Hq = 1.25 singular isothermal toroid. 



pole-on (6 = 90°), where p m i n reaches zero. Correspondingly, the isophotal aspect ratio 
r = b/a increases from a minimum value, that depends on the adopted model, to unity. 
For the three models considered in this paper, the minimum value of the aspect ratio is 
r m in ~ 0.6 for the Ho = 0.2 toroid, r m j n ~ 0.4 for the Ho = 0.5 toroid, and r m i n « 0.2 
for the Hq = 1.25 toroid. 

In Fig. [7| we show the maximum and minimum values of polarization as function 
of the aspect ratio r for the H = 1.25 toroid at A = 850 /im. The inverse correlation 
shown by the Figure between polarization degree and core aspect ratio is a characteristic 
signature of magnetic field configurations dominated by a poloidal component, and can be 
compared to observations to test the relative importance of poloidal vs. toroidal magnetic 
field components in dense clouds. In this sense, it is interesting to notice that of the three 
cloud cores represented in Fig. 1, the one showing the smallest range of variation (and 
absolute values) of p (Bl-c), is also the one characterized by the largest aspect ratio 
(r w 1, compared to r « 0.5 for L183 and L1544). A statistical analysis of the available 
observations may reveal the presence of a correlation between the range of polarization 
degree and the aspect ratio of cloud cores. As also shown by Fig. the maximum 



depolarization effect, measured by p„ 



Pn 



is obtained when r 



0.5 {9 



40°), and corresponds to a polarization reduction of a factor ~ 2, from p n 
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Fig. 6. Polarization degree as function of intensity at 850 /im (normalized to the peak 
value I max ) for the H = 0.2, 0.5 and 1.25 toroids for inclinations 6 = 20°, 40° and 60°. 
The multiple tracks visible in each panel are an artifact due to the sampling 
of the polarization degree over a regular square grid in the polarization maps. 



to Pmin ~ 7%. We notice that a detectable decrease in the polarization degree, say 
p maK — Pmin > 5%, occurs for inclinations ranging from ~ 20° to ~ 65°. For random 
orientations of the toroids, this interval corresponds to about 50% of all possible cases. 

3.3. Dependence of depolarization on wavelength 

In Fig. |H1 we compare polarization-intensity diagrams at A = 850 /im and A = 450 /im, 
obtained for the Hq = 1.25 toroid at 9 = 40°. As shown by the figure, the minimum 
and maximum value of polarization remain the same at these two wavelengths, but 
their dependence on intensity is different, the decrease of p with / being steeper at 
the longer wavelength. This effect can be easily understood: since the dust emission at 
A = 450 /xm is less concentrated than at A — 850 jum (see e.g. GGW), a given value of 
I/I max corresponds to a larger distance from the centre of the core (the intensity peak) 
at the shorter wavelength. Therefore the "polarization hole" appears restricted to higher 
values of I /Imax at A = 450 /xm. 
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Fig. 7. Variation of maximum (solid line) and minimum (dotted line) polarization degree 
as function of the isophotal aspect ratio r = b/a at A = 850 /im for the Ho — 1.25 toroid. 
The upper scale shows the inclination angle 9 (in degrees) of the line of sight with respect 
to the equatorial plane of the toroid. 
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Fig. 8. Polarization degree as function of intensity at 850 /.im (left panel) and 450 /im 
(right panel), normalized to the peak value J ma x, for the Hq = 1.25 toroid at 9 — 40°. 
Notice the different extent of the "polarization hole" , here identified with the region with 
low and uniform values of p on the right side of the dashed lines: within ~ 80 % and 
~ 50 % of the peak intensity at 850 and 450 /Ltm, respectively.) 



3.4. Effects of a non isothermal dust temperature distribution 

Fig. El shows the importance of computing the dust temperature distribution resulting 
from the heating of the external ISRF instead of assuming an uniform dust temperature 
through the cloud. In fact, the outwardly increase of the dust temperature from ~ 8 K 
to ~ 15 K for an externally heated cloud (see GGW) enhances, in the integration along 
the line of sight, the contribution to the Stokes parameters of the external layers of the 
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Fig. 9. Polarization degree as function of intensity at 850 /Ltm (normalized to the peak 
value /max) for the Hq = 1.25 toroid at 9 = 40°, in the case of external heating by 
the ISRF (left panel) and assuming an isothermal dust temperature distribution (right 
panel). 

cloud, where, as we have seen, the component of the magnetic field in the plane of the sky 
is generally larger than near the core's centre. Thus, an outward increasing temperature 
gradient contributes to the observability of the depolarization effect. In fact, given the 
current sensitivity of polarimeters at submillimetre wavelengths, that allow measurements 
of p only for / <; 0.2/ max (see Fig.^l, the value of p would be quite uniform through the 
sampled region if the cloud were isothermal (see Fig. El . Thus, the decrease of polarization 
shown by our model is actually the result of the combination of pinched magnetic field 
hourglass and a dust temperature gradient increasing outward. 

3.5. Effects of turbulence 

Finally, we remark that the predicted p—I diagrams of condensations formed in simu- 
lations of turbulent flows (Padoan et al. 2001) differ significantly from our Fig. as 
they show a large population of low-p and low-/ data points, which, in general, is not 
observed (cf. Fig.Q. Thus, at least in principle, p-I diagrams offer a way to distinguish 
between different explanations of the observed behavior of polarization in cloud cores, 
especially given the increased sensitivity expected from the next generation of detectors 
and polarimeters. 

In addition, our results show that a large-scale magnetic field with moderate equato- 
rial pinching can produce a significant deviation of polarization angles from the direction 
of the cloud's minor axis, up to about ±15° (see Fig. 3, 4 and 5). For example, in the 
case of the H = 0.2 toroid observed at an inclination 9 = 20° (Fig. 2), the distribution 
of polarization angles over the whole map is peaked on the direction of the cloud's minor 
axis (a consequence of the assumed axial symmetry of the model), but the dispersion in 
polarization angles around the mean direction is a x w 10°. This casts doubt on the use 
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of the Chandrasekhar- Fermi formula to estimate the magnetic field strength in molecular 
cloud cores, as this formula assumes that observed deviations of polarization angles from 
a given direction (of the order of a x w 10°— 15° in starless cores, see e.g. Crutcher et al. 
2004) are solely due to the presence of a turbulent (or better, "wavy" ) component of the 
field. 

4. Conclusions 

In this paper we have computed the polarization of the radiation emitted by dust grains in 
molecular cloud cores represented as magnetically supported equilibrium configurations. 
To this end, we have adopted the magnetostatic models of Li & Shu (1996) and the 
radiative transfer method developed by GGW. The Stokes parameters (and therefore 
the polarization degree) have been computed at two wavelengths (A = 850 ^im and 
A = 450 /um) and for various inclinations of the cloud with respect to the line of sight, 
assuming that the dust grains are elongated and aligned with the large-scale magnetic 
field. 

The main result of this paper is the demonstration that a significant depolarization 
effect, with characteristics very similar to those observed in actual cloud cores, can arise 
only because of geometrical effects, if the large-scale magnetic field has the equatorially 
pinched morphology predicted by magnetically dominated models. However, we do not 
claim that the field geometry is the sole means by which such an effect could be produced. 
There may also be contributions due to grain growth in the densest parts of a cloud, 
and turbulence. We also note that the assumption of axial symmetry of the density 
and field distribution implies that the model polarization vectors (after a 90° rotation), 
are symmetrically distributed with respect to the cloud's minor axis, at variance with 
observational evidence for some cores (see e.g. Ward-Thompson et al. 2000). A toroidal 
component of the magnetic field, ignored in the present analysis, can in principle account 
for the observed misalignment between the core apparent elongation and the average 
polarization position angle (see e.g. Vallee, Greaves & Fiege 2003). 
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